Gaussian processes may not be at the center of the current machine learning hype, but they are still used at the forefront of research. For instance, they were employed to automatically tune the MCTS hyperparameters for AlphaGo Zero. This type of nonparametric Bayesian algorithm can be very easy to use while providing rich modeling capacity and uncertainty estimates. However, the algorithm is not that easy to grasp, especially if you are new to nonparametric models. This section is proposed to present a more visual and intuitive introduction to Gaussian processes without totally abandoning the theory.
Introduction to Gaussian Process
A Gaussian Process (GP) is a powerful model that can be used to represent a "distribution (Gaussian) of functions" [65]. In a 2D regression problem, this distribution of functions can be visualized as a Gaussian distribution of curves, as shown in Fig. 7.5. That is, after the "training process" is finished, we would obtain an infinite number of fitting curves that are associated with different probabilities at a finite number of data points. Among them, we have a curve with the highest probability, i.e., the mean curve or the curve connecting all the mean values at the data points. We can also obtain confidence intervals, e.g., a range where a certain percentage of curves will fall in.
Figure 7.5: Conceptual illustration of Gaussian Process
Usually, we do the regression at a finite number of given data points (training and testing). In GP, the prediction or the label of each data point can be viewed as a random variable that obeys the Gaussian distribution. Accordingly, all the random variables corresponding to all the points form a multivariate Gaussian distribution, as illustrated in the following equation.
where ∼N( vec(mu), bar(Sigma))\sim \mathcal{N}(\vec{\mu}, \bar{\Sigma}) means following a Gaussian (normal) distribution characterized by the vector of means vec(mu)\vec{\mu} and a covariance matrix bar(Sigma)\bar{\Sigma}. Usually, the vector of means is assumed to be a vector of zeros, which is proven adequate to give the 'fitting function' or model enough flexibility. The overall location and trend of the curve are determined by the training data points.
The covariance matrix of the multivariate distribution is determined by a kernel function that defines the covariance coefficient between any two random variables in the matrix. The kernel function or the covariance coefficients calculated with it determine how the values between neighboring points are related to each other - the smoothness of the curve. A typical option for such a kernel function is a squared exponential function.
Gaussian process is different from many supervised learning methods, because it does not contain the typical training and testing steps, which is called eager learning. By contrast, it is similar to KNN , which is also a nonparametric algorithm, in that both of them belong to lazy learning. In such lazy learning algorithms, there is no obvious training step/stage, and most computing occurs during the testing/prediction. Because of this, this type of machine learning algorithm has a high demand for memory for storing the training data. This difference is very obvious in the implementation of such methods.
Using a 2D regression problem as an example, we have X_train, y_train, and X_test and need to predict y_test. The covariance matrix can then be obtained with these given data and the equations for multivariate Gaussian distributions. With the covariance matrix and a Zero vector of means, we can predict the possible curves consisting of the predicted values of the random variables at the testing data points. With a number of predicted curves, we can get the most possible (mean) curve and confidence intervals. Thus, there are training and testing data but no distinct steps for training and testing.
Also, as described above, the predictions were obtained based on the analytical equations from the training data (X_train, y_train) and X_test directly without assuming any specific math functions for the model. Instead, we only assume the values at all the data points follow a multivariate Gaussian function with a specific kernel function for the covariance coefficients. This is why Gaussian process is usually believed to be a nonparametric machine learning method.
Modeling Functions using Multivariate Gaussian
The key idea behind GP is that a mapping function can be modeled using an infinite dimensional multivariate Gaussian distribution. In other words, every data point in the input space is associated with a random variable, and the joint distribution of these variables is modeled as a multivariate Gaussian. This may not be easy to understand for people who are not familiar with multivariate Gaussian distributions. Let us use a simpler case, i.e., a unit 2D Gaussian, to see what it looks like.
First, let us take a look at a very simple multivariate Gaussian distribution that consists of multiple univariate standard Gaussian distributions. Such a function has a math formulation as follows:
Next, let us implement such a function. It is a probability distribution. So, we can sample 10 curves at 20 data points. Accordingly, II is 20 in this example.
import numpy as np
import matplotlib.pyplot as plt
def plot_unit_gaussian_samples(D):
p = plt.figure()
xs = np.linspace(0, 1, D)
for color in range(10): # 10 for 10 curves
ys = np.random.multivariate_normal(np.zeros(D), np.eye(D))
plt.plot(xs, ys)
return p
plot_unit_gaussian_samples(20)
As can be seen in Fig. 7.6, 10 curves sampled using this simple multivariate Gaussian distribution are illustrated using different colors. Each predicted value at any of the 20 data points is one sample from a univariate Gaussian distribution. One curve consisting of 20 random variable values was generated at one time. However, the curve is very 'noisy' or 'rough', which is much different from smooth math functions that we typically use for machine learning models. In other words, any two predicted values in each curve are not related. This is caused by the zero non-diagonal values in the covariance matrix, making the 20 predicted values in one curve not different from 20 predicted values from 20 separate standard Gaussian distributions.
Figure 7.6: Samples from a 20D Gaussian without kernel smoothing
In order to obtain smooth functions, we can generate a complicated covariance matrix using a squared exponential kernel equation. This kernel function, together with the Zero array for the mean vector, can be defined using the following code:
def kernel(X1,X2):
Sigma = np.empty((X1.shape[0],X2.shape[0]))
for i in range(X1.shape[0]):
for j in range(X2.shape[0]):
Sigma[i,j] = np.exp(-(X1[i]-X2[j])**2/2)
return Sigma
def mean(X):
Mu = np.zeros(X.shape[0]) #np.mean(X,axis=0)
return Mu
So, to get the smoothness we want, we will consider two random variables y_(i)y_{i} and y_(j)y_{j} plotted at x_(i)x_{i} and x_(j)x_{j} to have a covariance as cov(y_(i),y_(j))=k(x_(i),x_(j))\operatorname{cov}\left(y_{i}, y_{j}\right)=k\left(x_{i}, x_{j}\right). This dictates that the closer the two points, the higher their covariance.
Using the kernel function from above, we can get this matrix with "kernel(xs, xs)". Now, let us plot another 10 samples from the 20D Gaussian with the new covariance matrix. For this purpose, we will need to use mean(xs) and kernel(xs, xs) to replace np.zeros(D) and np.eye(D) in the above code, respectively. Then, plot the curves again, and we get Fig. 7.7.
Figure 7.7: Samples from a 20D Gaussian with kernel smoothing
Making Predictions using a Prior and Observations
Now, we get to the heart of GPs. As mentioned, though there are no explicit training and testing steps, we still need to predict the target values (y_test) at the testing data points (X_test). In GP, the mapping y=f(x)y=f(x) is represented by a probabilistic model p( vec(y)∣ vec(x))p(\vec{y} \mid \vec{x}) using a multivariate normal:
where bar(K)=k( vec(x), vec(x))\overline{\mathbf{K}}=k(\vec{x}, \vec{x}), in which the kernel matrix bar(K)\bar{K} is calculated with the kernel function introduced above, and m( vec(x))= vec(0)m(\vec{x})=\overrightarrow{0}, in which m( vec(x))m(\vec{x}) is the function for obtaining mapping vec(x)\vec{x} to its corresponding vector of means.
We have some training data with inputs vec(x)\vec{x}, and outputs vec(y)=f( vec(x))\vec{y}=f(\vec{x}). Now let us say we have some new points vec(x)_(**)\vec{x}_{*} where we want to predict vec(y)_(**)=f( vec(x)_(**))\vec{y}_{*}=f\left(\vec{x}_{*}\right). Accordingly, we assume all the data, including both the training and testing data, follows the same distribution. Therefore, we assume these two types of data from one multivariate Gaussian distribution, p(( vec(y)), vec(y)_(**)∣( vec(x)), vec(x)_(**))p\left(\vec{y}, \vec{y}_{*} \mid \vec{x}, \vec{x}_{*}\right), as follows.
where bar(K)=k( vec(x), vec(x)), bar(K)_(**)=k(( vec(x)), vec(x)_(**))\overline{\mathbf{K}}=k(\vec{x}, \vec{x}), \overline{\mathbf{K}}_{*}=k\left(\vec{x}, \vec{x}_{*}\right) and bar(K)_(****)=k( vec(x)_(**), vec(x)_(**))\overline{\mathbf{K}}_{* *}=k\left(\vec{x}_{*}, \vec{x}_{*}\right). As before, we stick with a zero mean. Thus, we can use what we know, vec(x)\vec{x} (X_train) and vec(x)_(**)\vec{x}_{*} (X_test), to calculate the covariance matrix.
Next, we will need to predict p( vec(y)_(**)∣ vec(x)_(**))p\left(\vec{y}_{*} \mid \vec{x}_{*}\right). This requires us to calculate vec(mu)_(**)\vec{\mu}_{*} and vec(Sigma)_(**)\vec{\Sigma}_{*}.
Now we have a posterior distribution over vec(y)_(**)\vec{y}_{*} using a prior distribution and some observations!
Example
In this subsection, we will use one example to illustrate the above process for implementing Gaussian process. Let us first assume we want to model data that follows a 5th-order polynomial function as given in Eq. 7.47:
Fig. 7.8 gives the true solution (curve) and training data (dots).
Figure 7.8: True solution and training data for Gaussian Process example
The mean vector and covariance matrix can be generated for this problem using the following code.
# Kernel
def kernel(X1,X2):
# Sigma = np.exp(-(X1.reshape(X1.shape[0],-1) - X2)**2/2)
Sigma = np.empty((X1.shape[0],X2.shape[0]))
for i in range(X1.shape[0]):
for j in range(X2.shape[0]):
Sigma[i,j] = np.exp(-(X1[i]-X2[j])**2/2)
return Sigma
def mean(X):
Mu = np.zeros(X.shape[0]) #np.mean(X,axis=0)
return Mu
Sigma = kernel(X_train,X_train)
Mu = mean(X_train)
The following is not a necessary step for problem-solving. We just want to see what we will get from multivariate Gaussian sampling without relating the training data to the testing data (Fig. 7.9).
# Sampling for showing the Gaussian samples
Sigma_sampling = kernel(X_true,X_true)
Mu_sampling = mean(X_true)
for n_sampling in range(10):
Y_sample = np.random.multivariate_normal(Mu_sampling,Sigma_sampling) #with smoothing (Sigma with
kernel)
plt.figure(2)
plt.plot(X_true,y_sample)
Y_sample = np.random.multivariate_normal(Mu_sampling,np.eye(Mu_sampling.shape[0])) #without
smoothing (Sigma without kernel)
plt.figure(0)
plt.plot(X_true,y_sample)
Figure 7.9: Multivariate Gaussian sampling without relating training data to testing data
The next step is the key step for conditioning multivariate Gaussian, which will obtain the coefficients for the testing data.
# Prediction
X_pred = X_true
K_whole = kernel(np.hstack((X_train,X_pred)),np.hstack((X_train,X_pred))) # Sigma matrix for the
combined Gaussian distribution of [X_train, X_pred]
K = K_whole[0:X_train.shape[0],0:X_train.shape[0]]
K_ast = K_whole[:X_train.shape[0],X_train.shape[0]:]
K_ast_ast = K_whole[X_train.shape[0]:,X_train.shape[0]:]
M = Mu
M_ast = mean (X_pred)
Mu_pred = M_ast + ( K_ast.T @ np.linalg.inv(K) @ (y_train - Mu).reshape((y_train.size,1)) ).reshape(
M_ast.size)
Sigma_pred = K__ast_ast - K_ast.T @ np.linalg.inv(K) @ K_ast
With the GP model, we can then sample ten curves and plot the result.
plt.figure(3)
for n_sampling_pred in range(10):
y_pred = np.random.multivariate_normal(Mu_pred,Sigma_pred)
plt.plot(X_pred,y_pred) # Plot n (number of iterations) predicted curves
Finally, we can obtain the average curve and confidence intervals by assessing the curves that we sample.
plt.plot(X_pred,Mu_pred) # Predict average prediction
# sigma = np.std(X_pred) # Incorrect
sigma = Sigma_pred.diagonal()**0.5 # Correct
plt.fill_between(X_pred,Mu_pred-2*sigma,Mu_pred+2*sigma,color='yellow',alpha=0.9) # Show the region
with mu +- sigma (68%) (2*sigma is 95%, 3*sigma is 99.7%)
The final results are shown in Fig. 7.10.
Summary
After going through the above theory and implementation, we can draw the following key points about Gaussian process.
Figure 7.10: Final results of Gaussian Process analysis
At each point for prediction, which can be viewed as the points for interpolation, the prediction at this point is treated as a random variable whose value can be drawn from a multivariate Gaussian distribution.
One set of predicted target values at all the prediction data points represents one possibility from the multivariate Gaussian distribution.
The average curve of all the predicted curves is the mean (most possible) prediction. The ranges of different confidence can also be obtained. Thus, GP does not give out one definite prediction, but instead, gives out different sets of predictions (corresponding to different models) that correspond to different probabilities.
Prediction at different points (including the training points) in a set of predictions (one curve) are related to each other via the kernel function in the covariance matrix of multivariate Gaussian. This kernel function, which can use different math functions, is given to define how to smoothen the curve; Usually, the closer the prediction point is to a training point, the smaller the difference in their target values.
The training points form a multivariate Gauss distribution with the prediction points. Thus, their values determine the general trends and means of predictions at the neighboring prediction points.
Besides, please note that the selection of the kernel is made by human experts while the determination of the parameters can be done automatically by minimizing a loss term. This is the realm of Gaussian process regression. Finally, the handling of noisy data also deserves close attention, especially when we cannot get perfect samples of the hidden function. In this case, we need to factor this uncertainty into the model for better generalization.
In summary, the locations (X_train) and target values (y_train) of the training points, the selection of the kernel function, and the locations of the prediction data points (X_pred) determine what the predictions look like. The predictions are different curves with different probabilities.
